function VARbs = AC_doProxySVARbootstrap_single(VAR,nboot,clevel)
% Code based on Gertler and Karadi (2015)
% Wild Bootstrap: Uses original estimates and draws (with replacement) the residuals.

jj=1;

% Removes a constant
res = detrend(VAR.res,'constant');

% Bootstrap simulations (indexed by jj)
while jj<nboot+1

    % T x 1 dimension array of -1 or 1 random numbers with prob 0.5-0.5 (Rademacher distribution)
    rr = 1-2*(rand(VAR.T,1)>0.5);

    % T x n randomly choose the sign of the time T shocks (all)
    resb  = (res.*(rr*ones(1,VAR.n)))';

    % Preallocate matrix of endogenous variables
    varsb = zeros(VAR.p+VAR.T,VAR.n);

    % Initial values
    varsb(1:VAR.p,:)=VAR.vars(1:VAR.p,:);

    % Use baseline estimates to generate endogenous variables, using newly simulated shocks
    for j=VAR.p+1:VAR.p+VAR.T
        lvars = (varsb(j-1:-1:j-VAR.p,:))';
        varsb(j,:) = lvars(:)'*VAR.bet(1:VAR.p*VAR.n,:) + VAR.bet(VAR.p*VAR.n+1:end,:) + resb(:,j-VAR.p)';
    end
 
    VARBS           = VAR;
    VARBS.vars      = varsb;

    % Resample also the proxies
    VARBS.proxies   = [VAR.m.*(rr(VAR.T-VAR.T_m-VAR.T_m_end+1:VAR.T-VAR.T_m_end,1)*ones(1,size(VAR.proxies,2)))];

    % Estimate again the VAR, and compute IRFs
    VARBS = AB_doProxySVAR_single(VARBS);
    
    % Save results into VARbs (output of function)
    for j=1:VAR.k;
        irs = VARBS.irs(:,:,j);
        VARbs.irs(:,jj,j) = irs(:);
    end

    jj=jj+1;   

end
    
for j=1:VAR.k
    VARbs.irsL(:,:,j)=reshape(quantile(VARbs.irs(:,:,j)',(1-clevel/100)/2),VAR.irhor,size(VARBS.irs,2));
    VARbs.irsH(:,:,j)=reshape(quantile(VARbs.irs(:,:,j)',1-(1-clevel/100)/2),VAR.irhor,size(VARBS.irs,2));
end